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ABSTRACT 

The spectra of several high-redshift (z > 6) quasars have shown indications for a 
Gunn-Peterson (GP) damping wing, suggesting a substantial mean neutral hydrogen 
fraction (xj^^ > 0.03) in the z w 6 intergalactic medium (IGM). However, previous 
analyses assumed that the IGM was uniformly ionized outside of the quasar's HII 
region. Here we relax this assumption and model patchy reionization scenarios for a 
range of IGM and quasar parameters. Compared to uniform reionization, patchy reion- 
ization imprints a different average damping wing profile with an associated sightline- 
to-sightlinc scatter. We quantify the impact of these differences on the inferred s}^^, 
by fitting the spectra of three quasars: SDSS J1148-I-5251 {z = 6.4189), J1030-^0524 
{z = 6.308), and J1623-h3112 {z = 6.247). We find that the best-fit values of x^^^ in 
the patchy models agree well with the uniform case. More importantly, we confirm that 
the observed spectra favor the presence of a GP damping wing, with peak likelihoods 
decreasing by factors of > few - 10 when the spectra are modeled without a damping 
wing. We also find that the Lya absorption spectra, by themselves, cannot distinguish 
the damping wing in a relatively neutral IGM from a damping wing in a highly ionized 
IGM, caused either by an isolated neutral patch, or by a damped Lya absorber (DLA). 
However, neutral patches in a highly ionized universe (a;^'^ < 10~^), and DLAs with 
the large required column densities (A^hi > few x 10^°cm~^) are both rare. As a result, 
when we include reasonable prior probabilities for the line of sight (LOS) to intercept 
either a neutral patch or a DLA at the required distance of 40 — 60 comoving Mpc 
away from the quasar, we find strong lower limits on the neutral fraction in the IGM, 
> 0.1 (at 95% confidence). This supports earlier claims that a substantial global 
fraction of hydrogen in the z w 6 IGM is in neutral form. 

Key words: cosmology: cosmic background radiation - dark ages, reionization - 
early Universe - diffuse radiation - large scale structure of Universe - quasars: general 
- emission lines 



1 INTRODUCTION 



of the population of Lya emitting galaxies. However, the 



The epoch of reionization, when the radiation from early 
generations of astrophysical objects ionized the intergalactic 
medium (IGM), offers a wealth of information about cos- 
mological structure formation and about physical processes 
in the early Universe. A variety of observations in the past 
few years have provided valuable insight into this epoch, 
through measurements of the cosmic microwave background 
(CMB) polarization anisotropies, quasar spectra, the 
kinetic Sunyaev-Zel'dovich (kSZ) signal and the evolution 
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In particular, fully black Gunn-Peterson (GP) troughs 
(IGunn fc Petersonlll965l ). have been found in a number of 
z > 6 Sloan Digital Sky Survey (SDSS) quasars, as well 
as in the qu asar ULAS J1120-I-06 41 recently discovered at 
z = 7.085 ("Mortl ock et~al]|201ll ) in the UKIDDS survey 



(jLawrence et al.ii2007l '). The presence of these troughs sug- 
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gests that these sources may be re vealins the tail-end of 
rcionization (see, e.g., the review bv lFan. Carilli fc Keating 
[2OO& I. However, a neutral fraction as low cLS ^fjj 10"* 
in a patch of the IGM is sufficient to produce a black GP 
trough - i.e. to make the flux undetectable at the wave- 
lengths corresponding to resonant Lya absorption. As a re- 
sult, constraints inferred about global properties of reioniza- 
tion - such as the volume-averaged neutral hydrogen frac- 
tion - require modeling the spatial distribution and 
time-evolution of neutral hydrogen in the IGM (see e.g., 
lMesingej[2oTol ). 

A unique signature of significant neutral hydrogen in 
the IGM is the presence of a long tail of absorption, extend- 
ing to wavelengths far away from t he resonant GP troug h 
- the so-called GP damping wing (iMiralda-Escudel 1 19981 ) . 
The optical depth in the damping wing is 5 — 6 orders of 
magnitude lower than in the resonant core of the Lya line. 
Therefore this damping wing is detectable only if the neutral 
fraction is significant (x]§^ > few xlO"^). 

In the SDSS quasars with a fully black GP trough, 
flux remains detectable over a wavelength range of a few 
X 10 A blueward of the Lya line center, attributable to the 
nearly fully-ionized region surrounding the quasars. This re- 
gion is highly ionized, but still contains sufficient residual 
HI to produce observable resonant line-of-sight (LOS) ab- 
sorption, i.e. the Lyman a forest. Additionally, as the pho- 
tons pass through the (possibly) neutral IGM in the fore- 
ground, the absorption in the GP damping wing causes 
a further depression of the flux. These two sources of 
Lya absorption differ markedly: the resonant absorption 
within the ionized bubble strongly fluctuates with wave- 
length, whereas the absorption from the GP damping wing 
smoothly and monotonically increases blueward, with only 
a mild wavelength-dependence. This allows, in principle, 
the two different sources of opacity to be s eparated, given 
sufficiently high-quality quasar spectra (ICen fc HaimanI 
I2OOOI : iMadau fc Re3l2000l : iMesinger. Haiman fc Cenll2004l ). 
However, any robust interpretation requires modeling due 
to the near- degeneracy and uncertainty of the quasar 
and IGM properties (e.g., iMesinger fc Hainianl '2004'. '2007*; 
Bolton fc Haehneltl l2007 ej|bl: I Mesinger fc Furlanettq i2008i : 



MaseUi. Ferrara fc Galle rani 200^. 



In the case of the UKIDDS quasar ULAS J1120-H0641, 
the GP trough is located at a wavelength much closer to 
the Lya line center, i.e. the ionized region appears markedly 
smaller. As a result, the situation for this object is somewhat 
different from the SDSS quasars: the red GP damping wing 
can cause sig nificant (1 0-20%) absorption on the red side of 
the Lya line (|Bohon et a l. 2011)F1 Therefore, the GP damp- 
ing wing and the intrinsic quasar spectrum must be modeled 
simultaneously, and care must be taken to include the uncer- 
tainties in the latter (by comparison, the uncertainty in the 
intrinsic spectrum plays a sub-dominant role when deriving 
lower li mits on the neutral fracti on from the SDSS quasar 
spectra: iKramer fc HaimanI I2OO9I ). We defer the analysis of 
ULAS J1120-I-0641 to future work. 



In our previous an alysis of the SDSS quasars 
l|Mesinger fc HaimanI l2007l : hereafter MH07), we modeled 
the fluctuating resonant absorption within the ionized zone, 
detecting signatures of GP damping wing absorption in two 
out of three SDSS quasars: J1030-I-0524 and J1623-I-3112. 
Treating each as a free parameter, we were able to simul- 
taneously place limits on the mean neutral fraction Xj^^, 
the quasars' ionizing luminosities, and the sizes of the sur- 
rounding ionized regions. Most intriguingly, we statistically 
inferred the presence of a smooth absorption component cor- 
responding to the GP damping wing, setting strong lower 
limits x^^j^ > 0.03 for both J1030-h0524 and J1623-h3112 
(with best-fit values of Ihi^ = 1 in both cases). 

In MH07, a uniform ionizing background was assumed. 
However, the topology of reionization is likely very inho- 
mogeneous, with ioni zed bubbles growing around highly 
clust e red galaxies (e.g . Furlanetto. Zaldarriaga fc Hernquist 



|2004 iMcbuinn et al.! [20071 : iTrac fc CenI l2007l : IZahn et al 
2011b|). The shape of the GP damping wing in a universe 
with such a "Swiss-cheese" reionization topology differs 
from that in a smoothly ionized IGM (as could result from 
more exotic r eionization scenarios by X-ray photons; e.g. 
lHaiinanll201ll) 



1 The 2 = 6.44 quasar CFHQS J0210-0456 (|WiIIott et aI.ll2010a^ 
has a similarly small ionized region (f^l.7 Mpc) and the same 
comment applies to this source; however, this source is intrinsi- 
cally fainter and has a lower S/N spectrum. 



iMesinger fc Furlanettol (|2008l ) compared GP damping 
wing profiles in the se two different scenarios (see also 
iMcQuinn et al.ll2008l ). Their results suggest that the differ- 
ence could be significant, and that it can introduce a bias 
and scatter in neutral fraction constraints. In this paper, 
we use physical models for patchy reionization based o n 
a semi-numerical simulation (jMesinger fc Furlanettdl2007l ). 
and infer constraints on the neutral fraction and on quasar 
parameters, using the observed spectra of the same three 
quasars that were previously analyzed in MH07. Our main 
goal is to assess whether allowing for patchy reionization 
significantly modifies previous results. Apart from relaxing 
the assumption of a uniform ionizing background, our anal- 
ysis improves on MH07 by (i) using updated quasar redshift 
measureme nts and constraints o n the quasars' ionizing lu- 
minosities ijCalverlev et al.l201lf ). (ii) incorporating the fluc- 
tuating ionizing background from galaxies near the quasar, 
and (iii) a better statistical characterization of the confi- 
dence levels on the model parameters, explicitly contrasting 
results from a Bayesian analysis with a parametric boot- 
strapping procedure (as opposed to inferring probabilities 
from a Kolmogorov-Smirnov test as in MH07). 

The rest of this paper is organized as follows. In § [2j we 
describe our analysis technique, including a description of 
the patchy reionization simulations, as well as a brief sum- 
mary of how we generate mock absorption spectra and fit 
them to the observational data. In § [31 we present our re- 
sults. In §131 we discuss various aspects of our results, such as 
the origin of the constraints we obtain, and the robustness of 
our results, including the possible presence of damped Lya 
systems (DLAs) . In § (5] we summarize our key findings and 
present our conclusions. Throughout this paper, we adopt 
the background cosmological parameters (Oa, Om, i^b, n, 
0-8, Ho) = (0.73, 0.27, 0.0455, 1, 0.76, 70 km s"^ Mpc"^), 
cons istent with the seven-year results by the WMAP satel- 
lite l|Komatsdl2010l ). Unless stated otherwise, all quantities 
are quoted in comoving units. 
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2 ANALYSIS 

Our analysis has four distinct components: (i) running semi- 
numerical simulations to model patchy reionization, as well 
as the spatial distribution and mass function of halos, and 
the IGM density and velocity fields; (ii) using the simula- 
tions to create mock absorption spectra; (iii) modeling the 
quasars' intrinsic emission spectrum; and (iv) comparing the 
observed and simulated absorption statistics. We first dis- 
cuss our patchy reionization models (§ I2.1|l . Then in § 12.21 
we present our procedure for creating mock Lya absorption 
spectra (ii). Finally, we discuss the comparison of mock and 
observed spectra (iii-iv) in § 12.41 



2.1 Modeling Patchy Reionization 

We use the publicly- available , semi-numerical code Dexivjfl 
l|Mesinger fc Furlanettol boOTi ) to generate evolved density, 
halo, and ionization fields at z = 6.5. Here we briefly 
outline the p rocedure, and refer the interested reader to 
iMesinger fc Furlanettg (|2007f ) for more details. 

Our simulation box is L = 250 Mpc on a side, with 
a particle grid cell size of 0.14 Mpc. Halos are filtered out 
of the 1800^ linear density field using excursion-set theory. 
Halo locations are then mapped to Eulerian coordinates at 
a given redshift using first-order Lagrangian perturbation 
theory l|Zerdovichlll970h . The resulting halo fields have been 
shown to match both the mass function and statistical clus- 
tering properties of halos in N-body simulation s, well past 
the linear regime (jMesinger fc Furlanetto|[2007f ). 

The evolved (non-linear) density field is computed in 
the same manner, by perturbing the 1800^ Lagrangian den- 
sity field. The resulting particle locations are then binned 
onto an Eulerian 600'^ grid, and the corresponding evolved 
velocity field is re-computed from the density field. The fi- 
nal density and ionization fields thus have a grid cell size of 
0.42 Mpc, approximately corresponding to the Jeans length 
in the ionized IGM at mean density. The statistical prop- 
erties of the density and velocity fields have been shown 
to match those from a hydrodynamic simul ation remark- 
ably well (|Mesinger. Furlanetto fc CenI 120111 ) . Specifically, 
the agreement between the density fields at z « 7 is better 
than <10% for both: (i) the PDF (smoothed on comparable 
scales) for roughly a dex around the mean (accounting for 
the vast majority of the volume), and (ii) the power spectra 
at fe < few X Mpc~^. 

To generat e the ionization field, we use the excursion - 
set formalism (|Furlanetto. Zaldarriaga fc Hernguisd l2004h . 
which compares the number of ionizing photons produced 
in a region of a given scale to the number of neutral hy- 
drogen atoms inside that region. Specifically, we identify 
ionized cells in our box as those which meet the criterion 
/coll ^ C~^i where is an ionizing efficiency parameter and 
/coll is the collapse fraction smoothed around a cell at (x, z) 
on decreasing scales, R. The collapse fraction is computed 
using the resolved halo field, including halos down to a min- 
imum mass of 1.7 x 10* M©, corresponding to the atomic 
cooling threshold at 2 ~ 6. The resulting ionization mor- 
phologies match those generated with cosmological radiative 



http: / /homepage. sns.it/~mesinger/Sim. html 



transfer algorithms relatively well, with the power spectra 
agreeing to within ~ 10-20% down to the Nyquist frequency 
(.Zahn et al. 2oTTbh . 

We follow the procedure described in 
IMesinger fc Diikstral (|2008l ) to model the galactic compo- 
nent of the ionization rate, Pbg, inside the HII regions, 
assuming that each galaxies ionizing luminosity is propor- 
tional to its halo mass. At every cell location (x, 2;), we sum 
the contributions of halos with mass A'h at x^ : 



(l + zY 

rbg(x,2) = ; (TH^iou ^ 



Ah 



Itt 



I / 



-|x-Xi|/A 



(1) 



where gh is the Lyman limit cross-section, and A is the ioniz- 
ing photon mean free path inside HII regions. For simplicity 
and computational convenience, our fiducial analysis uses a 
constant mean free path, A = 60 (c omoving) Mpc, consisten t 
at 1-(T with observations at 2 ~ 6 jSongaila fc Cowi3l2010l ): 
we return to this point further below. The rate per unit mass 
at which ionizing photons are released into the IGM, eion, 
governs the overall normalization of eq. ([T l. This normal 



ization is chosen to match observations (see lCalverlev et al.l 
l201ll ). resulting in an average ionization rate per hydrogen 
atom of Fbg = 10~^^ s~^. In ot her words, we fix the mean 
galactic ionization rate to match [Calverlev et al.l (|201lh . but 
use individual halo locations to model the local spatial fluc- 
tuation^ around this mean, according to eq. ([TJ. 

We vary the ionizing efficiency parameter, i^, to 
generate a suite of ionization fields corresponding to 
various values of at 2 = 6. Specifically, we 

use the following eight different values: iin^ = 
{0.01, 0.06, 0.10, 0.38, 0.53, 0.64, 0.73, 0.88}. We then extract 
2 X 10'* lines of sight (LOSs) centered on the two most mas- 
sive halos in the simulation volume (Af ~ 3 x lO^^M©). The 
LOSs are drawn from randomly chosen directions. The re- 
sulting density, velocity and ionization LOSs are then used 
to construct mock Lya spectra, as described below. 



2.2 Modeling Quasar Absorption Spectra 

We compute the residual, local neutral fraction along each 
LOS in the HII region surrounding the quasar, whose extent 
we treat as a free parameter, Rs- Inside the ionized region, 
we assume ionization equilibrium with radiation from the 
quasar (Fq) and the background galaxies (Fbg), 



(Fbg + FQ)a;HinH = aBXcff?^H(l — a^Hi) 



(2) 



where nH(x,2) is the local (non-linear) number density of 
hydrogen atoms, xm is the local neutral hydrogen fraction, 
Xcff = 1-08 accounts for additional electrons from singly- 
ionized Helium (assuming that the Helium and Hydrogen 



For simplicity and to avoid introducing uncertain scalings we 
use the same UVB for each - i.e. we ignore the "shadow- 

ing" from the neutral islands in patchy reionization models with 
> 0. This can become important for large aij^'^, when the 
sizes of many ionized bubbles are smaller than the assumed mean 
free path of A = 60 Mpc, increasing the spatial fluctuations in 
the UVB. However, this is unlikely to influence our results, as the 
absorption statistics are much more sensitive to th e fluctuations 
in the density field jMesineer fc Furlanettol l2009al : see also i|4.5l 
and the associated discussion). 
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ionized fractions are equal with a Helium mass fraction of 
Yua = 0.24), and as = 1.7 x 10"^^cm^s"^ is the case-B re- 
combination coefficient. The latter is evaluated at the tem- 
peratu re T = 1.8 x 10^ K (consistent with recent measure- 
ments; iBolton et aLll2012l n. 

Equation [5] is solved at every step along the LOS, using 
the non-linear density from our simulation. The contribu- 
tion to the ionization rate from the quasar drops off approx- 
imately as the inverse square of the distance, and is given 
by 



To = 



1.6 X W^\fr 



dv 



(-)■ 



l + Z 

1 + 2Q 



Here ct^ ~ 6.3 x 10~^*cm'^(z//!/H) is the frequency- 
dependent hydrogen ionization cross-section, h is Planck's 
constant, zq is the observed redshift of the quasar, z is the 
total redshift corresponding to a pixel along the LOS, includ- 
ing pecuhar velocity [taking (1 -I- 2) = (1 -I- 2:o)(l + ""pec/c), 
where zo is the redshift that would correspond to the Hub- 
ble flow at the comoving distance of the pixel], = 
3.29 X lO^^Hz is the ionization threshold of hydrogen, /r 
is the quasar's ionizing luminosity in units of 1.3 x lO^'^ s~^ 
(following MH07), and dh = dL{zo,zQ) is the luminosity 
dist ance between the quasar and the for eground pixel (see, 
e.g.. IPhillipps. Horleston fc Whitd[20o3 ). Equation © fur- 
ther assumes an intrinsic power-law spectrum for the 
ionizing radiation of the qu asar, matching standard quasar 
templates (e.g. iTelfer et al.ll2002l '). 

As in MH07, our model has three free parame- 
ters: (1) the LOS distance to the edge of the HII 
region surrounding the quasar, Rs, (2) the quasar's 
ionizing luminosity parametrized by /r, and (3) the 
mean volume-weighted neutral fraction in the IGM, 
x^'^. Our model grid covers a different range of 
Rs for each quasar: in units of comoving Mpc, for 
J1148-h5251, Rs = {37,39,41,43,45,47}, for J1030-I-0524, 
Rs = {52,54,56,58,60,62}, and for J1623+3112, Rs = 
{37,39,41,43,45,47}. In each case, we cover the range 
fr = {0.1,0.4,0.7,1.0,1.3,1.6,1.9,4.0,6.0,8.0,10.0}, and 
the eight different values of listed above. 

It is worth pointing out that, in general, for a fixed 
neutral fraction, Rs could be determined by the ionizing 
emissivity and the lifetime of the quasar, the surrounding 
galaxies, or a combination of both. As such, Rs may be 
an integral , indirect measure of galaxy and quasar param- 
eters (e.g.. iWvithe fc Loebl [20071 ) . but here we simply take 
it to directly describe the ionization topology surrounding 
the quasar. We chose a range of Rs values greater than or 
equal to the comoving distance between the redshift of the 

* While case-B recombination, corresponding to optically thick 
media, is typically assumed in the reionization literature, the ap- 
propriate coefficient is somewhat uncertain. Depending on the 
spectrum of the ionizing sources and the abundance and density 
profiles of optically thick Lyman limit systems, the case-A value 
could be more appropriate (see discussion in iMiralda-Escudj 
120031) . Instead of attempting to model the appropriate coeffi- 
cient, we note that its value is degenerate with both the assumed 
temperature, and also with the ionizing luminosity. In particu- 
lar, the diflFerence between the two coefficients is a factor of ^ 1.6 
llOsterbrocklll974l) , well within the range of possible temperatures 
in quasar HII regions (Bolton et al. 2012). 



quasar and the redshift of the first (i.e. reddest) pixel in the 
quasar's GP-trough. As we utili ze the recently publi shed 
revised redshift measurements of ICalverlev et al]|201ll . our 
Rs values are somewhat different (i.e., larger) than those in 
MH07. The wavelengths corresponding to our Rs grid are 
indicated by the (red) vertical dashed lines in Figure [T] 

For each given combination of (_Rs, /r, Spji^), we cre- 
ate a mock absorption spectrum corresponding to each LOS 
extracted from our simulation box. The spectrum is sim- 
ply assumed to be black at wavelengths corresponding to 
R > Rs and is effectively ignored from our analysis □ At 
(•3^ R < Rs, the total optical depth includes resonant absorp- 
tion from neutral HI inside the ionized region, as well as 
absorption through the damping wing by neutral hydrogen 
located outside R > Rs- We first create an array of ob- 
served wavelength bins Aj, and for each bin, we sum the 
optical depth from every pixel i along the LOS, 



nieC 



Ss 



1 + zo 



+ Zi 



(4) 



where e = 4.8 x 10"^" esu is the charge of the electron, 
fa = 0.4162 is the Lya oscillator strength, rrie = 9.1 x 10~^* 
g is the mass of the electron, Zi is the total redshift of pixel 
i, including its peculiar velocity, Xi is the neutral fraction 
and rii is the total hydrogen density in pixel i, Ss/{1 + zo,i) 
is the (physical, as opposed to comoving) width of the pixel, 

and 0a is the normalized absorption profile. 

The frequency-dependence of the Lya line is (jPeeblesI 

1 19931 . eq. 23.97) 



0H = 



(5) 



where A^ = 6.25 x 10*s~^ is the radiative decay rate. A com- 
mon appro ach in the liter ature is to use the convolution of a 
Lorentzian ()PeebleJl993l . eq. 23.43) with the thermal broad - 
enmg, i.e. a Voigt profile ()Press. Rvbicki fc Schneide3ll993h 
near the resonance oj « Wq and an analytic solution for the 
damping wing far away from resonance that ignores Doppler 
broa dening but includes the {u)/uJaY Rayleigh-scattering 
tail (|Miralda-Escudel Il998l ) . As our results may depend on 
intermediate regimes between the damping wing and line 
center, in principle, the thermally broadened line profile 
should be computed by convolving equation ([5} with the 
temperature-dependent Doppler profile. In practice, we have 
approximated this full convolution with 



0a ^ 0voigt 



(6) 



where 0voigt is the Voigt profile for our tem perature T . This 
appro ximation approaches that given by iMiralda-Escudel 
(199a) on the far red side of the damping wing. 



^ In principle, spikes of transmission can appear at > Ks-, 
outside the quasar's zone of influence. The statistics of the 
spectrum in this range can also place constraints on the ion- 
izing backgr ound and on the volume-filling fraction of neu- 
tral patches llCroftlll9 98': 'Songaila fc Cowiell20ol lBarkana|[20o3: 
iGallerani. Choudhurv fc Ferrara 200(|), although in practice, con- 
straints on the latter require an independent measurement of the 
UVB (Mesinger 2010). We emphasize that our results below are 
independent of these statistics and arise solely from the GP damp- 
ing wing. 
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In our patchy reionization models, specifying de- 
fines the ionization morphology. It is unlikely that the edge 
of the surrounding HII region will happen to lie at any given 
fixed distance, Rs- Therefore, in order to keep Rs as an in- 
dependent free parameter, we "shift" the ionization array 
(keeping the density array along the LOS fixed) by the min- 
imal amount so that a neutral pixel falls at the location 
at Rs- However, in the highly-ionized models, we often en- 
counter LOSs without any neutral pixels. In these cases, we 
keep the density and velocity field along the LOS, but we 
"recycle" the ionization array from another, randomly cho- 
sen LOS to make sure that a neutral pixel is located at Rs- 
This shifting is done in order to treat Rs as a free param- 
eter, while preserving the ionization morphology from the 
simulation at i? > Rs (which determines the GP damping 
wing). After presenting these conservative constraints, we 
then combine them with the (model-dependent) probability 
that a neutral pixel is, in fact, found at the radius Rs (see 
the discussion in i]3.2l below). This takes into account that 
a more ionized universe generally has larger HII regions by 
introducing a prior on Rs- 



2.3 Model Parameters 

In summary, our fiducial analysis has three free parameters 
(analogous to MH07) ; 

• Xui^ ~ the mean volume-weighted neutral fraction in 
the IGM. At a given redshift, is a function of ^ and 
Tvir (see belowfl. 

• /r - the quasar's ionizing luminosity in units of 
1.3 X 10^'' s~^ 

• Rs - the LOS distance to the edge of the HII region 
surrounding the quasar. 

Our models include additional parameters, which we either 
do not vary or which are used to derive the above quantities. 
Ideally one would like to explore as large a parameter space 
as possible; however, in practice, computational costs have 
limited us to vary at most three parameters simultaneously. 
In Sections 4.4-4.6 below, we perform additional investiga- 
tions to assess the robustness of our results to degeneracies 
missing in our three-parameter models. The additional pa- 
rameters include: 

• (" - the ionizing efficiency of high-redshift galaxies. 
This quantity can be defined as ( — fcsc,f*N.y /{I + riroc), 
where /esc is the fraction of ionizing photons produced by 
stars that escape into the intergalactic medium (IGM), /* 
is the star formation efficiency, is the number of ionizing 
photons per stellar baryon, and rircc is the mean number 
of recombinations per baryon. For reference, fasc ~ 0.1, 
/* — 0.1, A''.^ = 4000 (appropriate for PopII stars), and 
riioc ~ 1 yield — 20. As described above, we vary ^ to 
generate ionization maps, computing the resulting xj^^. 



^ In principle, the ionizing morphology at fixed can also be 

a functio n of A; however, this d ependence is very weak for A > 10 
Mpc (e.g. lMcQuinn et~al]|2007l ') . ~ 



• 2vir - the minimum virial temperature of halos hosting 
star-forming galaxies. As is common in reionization litera- 
ture, we set this to the atomic-cooling threshold, Tvir = 10* 
K. Although Tvir could have a somewhat different value, 
the ionization morphology at fixed (the quantity 
of interest in this work) is fairly r obust to (reasonable , 
physically-motivated) changes in Tyi r (jMcQuinn et al.ll2007l : 
iMesinger. McOuinn fc SpergelllioTll ). 

• fbg - the mean value of the galactic UVB. As men- 
tioned above, we set this value to fbg — 10^^'^ s~^. 
Although this m atches observations o f the a ~ 6 Lyman 
alpha forest (e.g. ICalverlev et al.ll201ll ). such measurements 
have large uncertainties. Luckily, our conclusions are not 
sensitive to the UVB, as we discuss below. 

• A - the ionizing photon mean free path inside HII re- 
gions, generally set by Lyman limit systems (LLSs). We use 
a fiducial value of A = 60 (como ving) Mpc, consistent a t 
1-0" with observations at 2 ~ 6 (jSongaila fc Cowi3 12OI0I ). 
The mean free path is used to compute the fluctuations 
in Fbg, the galactic UVlfl. In principle, A is poorly con- 
strai ned at high redshif ts, and could even vary spatially 
fe.g. ICrociani et al.ll201lh . However, our conclusions are not 
very sensitive to Fbg, for several reasons. First, for most 
of the relevant spectral range and parameter choices, the 
quasar's flux dominates the ionizing background, with galax- 
ies becoming important only in the relatively unimportant 
high-_Rs, low-/r, low-a;J^^ regime ( see more discussio n of 
this point below in ij | 4.5I and also iLidz et al ] l2007bl and 
iBolton fc Haehneltll2007al ). Second, the profile of Fbg (as a 
function of the distance from the quasar) is not degenerate 
with any of our three free parameter^. Finally, as mentioned 
above, the average transmission is far more sensitive to the 
density fi eld than to inhomogen eities in the background fiux 
(Mesinger fc Furlanettoll2009al ). 



2.4 Fitting the free parameters in our model 

2.4-1 Pixel Optical Depth Distributions 

In order to compare our simulated spectra to the observa- 
tions, we also need an intrinsic (unabsorbed) quasar emis- 
sion template. These are obtain ed by fitting the K eck ESI 
observations of the quasars (e.g.. lWhite et a h 2003) only on 
the red side of the Lya emission line, where IGM absorption 
is minimal. The intrinsic spectral fits are taken from MH07, 
to which we refer the interested reader for details. Briefly, 

^ It is also possible that LLSs along the line of sight from the 
QSO attenuate the quasar's flux, Fq so that it deviates from 
the canonical oc r~-^ behaviour. Such LLSs would need to have 
column densities low enough to escape detection in HIRES and 
ESI. We return to this issue in il4.6l 

* Fbg is spatially fluctuating around the mean value, with 
an expected va riance on the sc a le of our cell size of a 
factor of ~few llLidz et al.l )2007al : iMesineer fc Diikstral I2OO8I: 
iMesinger fc Furlanettdl2009bl )" There is a slight evolution of the 



mean with distance from the quasar, due to the biased nature of 
the clustered galaxies (which is accounted for in our models). The 
fact that none of our free parameters have similar imprints on the 
QSO spectra suggests that any inaccuracies in our modeling of 
Fbg are unlikely to have a strong bias on our results. 
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the model fit to the data consists of the sum of a power-law 
continuum, a double-Gaussian Lya emission line, and a NV 
emission line. The observed continuum-normalized flux in 
wavelength bin j is then obtained as 

^obs,j ~ Fobs,j / Fcont^j , (7) 

where Fcontj and Fobs,], are the intrinsic and observed flux, 
respectively. We compare this observational signal with our 
mock continuum-normalized spectra, given by 

Fsim,] = exp(~rj). (8) 

Simulated optical depths can be far in excess of what is ob- 
servable, so to ensure a fair comparison between simulation 
and observation, we impose a floor of minimu m flux given 
by th e observational errors {jJoba.j, taken from IWhite et all 
I2OO3I ): 

FminJ = ^'^obs,j/F^^^^^^, (9) 

with the constant chosen to be y4 = 3 for our analysis, cor- 
responding to 3-cr detections. For all values of Fj < Fmin.j 
we treat the simulated and observed continuum-normalized 
flux to be the same as zero flux. This floor corresponds 
to a maximum observable optical depth ranging between 
Tob3 = 4.75 to 5.5. The continuum-normalized spectra of 
the three quasars, subject to these conditions, are repro- 
duced in Figure[T] Since present-day simulations do not have 
the dynamic range to model biases surrounding the bright, 
rare quasars, while at the same time resolving the Lyman- 
Q forest, we excise the region 25 A blue- ward of the Lya 
peak from our analysis. This roughly corresponds to the ex- 
pected mean radius (~ 10 comoving Mpc) of the large-scale 
overd ensity surrounding such quasars (e.g. iBarkana fc Loebl 
|2004| ). 

We then compare the observed spectra to simulated 
mock spectra, using the tw o-sided Kolmogo rov-Smirnov 
(KS) distance statistic Dks (|Press et al.lll993 ). applied to 
the observed vs. simulated flux density PDFs. Each spec- 
trum was binned into three wide wavelength ranges, with 
two of the bins chosen to lie red- ward of the pixels where the 
observed spectrum is consistent with zero flux - nominally 
corresponding to the apparent edge of the ionized bubble 
surrounding the quasar. The apparent edge falls within the 
third, bluest bin (see Figure [T]). Within each of these three 
bins, the PDFs of the optical depth was obtained, for a given 
combination of {Rs, /r, 2;hi^)i using all LOSs and all pixels 
within the bin. Each of the three model PDFs were then 
compared to the corresponding optical depth PDFs inferred 
from the observed spectra. 



To address the impact of all of these issues, we explic- 
itly compute the PDF of the KS distance statistic Dks in 
each of our models. We flrst compute the pixel optical depth 
PDF for the entire set of 2 x lO'' LOSs - this constitutes our 
best estimate for the t PDF predicted in a given model. We 
then compare this prediction to the r PDF for each individ- 
ual LOSs in the same model, thus obtaining the (cumula- 
tive) PDF of Dks in this model. We have found that these 
CPDFs diff e r sign iflcantly from the fitting formula Qks in 
iPress et al.l (|l992l '). We traced these differences to the two 
reasons (i) and (ii) above. First, the fioor imposed by equa- 
tion ([5} causes a pile-up of nearly identical optical depths 
(corresponding to the number of "black" pixels) in the r 
PDF. This produces sharp steps in the sample CPDFs, and 
this unusual feature modifies the PDF of Dks (this effect is 
similar to a decrease in the number of effective data points). 
Second, correlations between neighboring pixels skew the 
CPDF away from Qks- We demonstrated this by first mask- 
ing out the black pixels from the analysis, and then picking 
r values from random LOSs to generate 2 x 10* uncorrelated 
r distributions. In this case, the CPDF of Dks indeed was 
consistent with Qks{Dks)- 

For each {Rs, fr, im^) model (and in each of the three 
spectral bins i of each quasar) we measure the KS distance 
DKS,obs, by comparing the r CPDFs in the given model 
(averaged over all LOSs) with the corresponding r CPDF 
in the quasar spectrum. We then construct the expected 
CPDF, P/l" \ -igmX> Dks), of the KS distance in this 

v-^s ' jr ''^Hi ' 

model, by comparing the r CPDF of each LOS to the mean 
r CPDF constructed from all LOSs. The probability that 
Dks exceeds DKS,ob3, ^(fisjj,,siGM)(> DKS,obs), is assigned 
to this model. Finally, multiplying the three probabilities in 
the three wavelength bins gives the model's overall likeli- 
hood p{Rs,fr,xW,^) = nLi^(«;.Vr,^jGM)(> DKS,obs) of 
the model. The best-fit model was identified as the one that 
maximizes this overall likelihood. 

We next place confidence limits around this best-fit 
model by interpreting p(J?s, /r, Sjn^) directly as Bayesian 
probability densities, and integrating these probability den- 
sities over our parameter space to obtain the confidence in- 
tervals. The results from this approach depend on the (ar- 
bitrary) parameterization of the parameter space volume. 
Here we chose Rs , log fr and log - effectively adopt- 

ing a uniform prior on these quantities. As a check on the 
robustness of these Bayesian confidence limits, we compare 
these results to those inferred from a parametric bootstrap- 
ping procedure (see § 14. 2[) . 



2.4-2 Best-Fit Models and Confidence Limits 

In MH07, for each set of model parameters, the product of 
the three KS probabilities associated with the three Dks 
values, was used as the likelihood of that model. The prob- 
ability Qks{Dks) is given by the commonly us ed approx- 
imate analytic fitting formula jPress et al.lll993 ). However, 
this approach is potentially problematic, since: (i) the an- 
alytic formula is only accurate for "well-behaved" distribu- 
tions; (ii) r values in nearby pixels can be correlated, and 
(iii) the r values in each of the ~10 consecutive pixels within 
a bin are drawn from slightly different underlying probabil- 
ity distributions. 



3 RESULTS 

3.1 Uniform Reionization 

Before attempting to improve on the results in MH07, we 
first reproduced their main results, assuming, as in that pa- 
per, that the IGM outside Rs is ionized by a spatially uni- 
form background. A few of the details of our analysis, how- 
ever, still differ from MH07: 

(1) To accommodate the revisio ns in the redshifts o f the 
three quasars (by up to « 0.03; ICarilli et al.ll2O10l ). we 
modify the boundaries of our three wavelength bins com- 
pared to those used by MH07. Likewise, we slightly adjust 
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and expand the set of values for Rs used in our grid of mod- 
els. In particular, for J1030+0524, which has a Ly/3 trough 
that lies at a significantly larger distance from the quasar 
than the Lya trough, we impose a minimum size for the ion- 
ized region to i? > 52 Mpc, corresponding to the location of 
the Ly/3 trough (as opposed to 40 Mpc in MH07, using the 
old redshift). 

(2) Rather than having a uniform UVB inside the 
quasar's HII region {R < -Rs), we compute the inhomo- 
geneous UVB from the nearby galaxies. 

(3) We find the best-fit models using our own deter- 
mination of the KS probabilities, rather than the standard 
fitting formula, as discussed in § 12.4.21 above. 

(4) We add a strong prior for the probability distribu- 
tion of the quasar's ionizing flux / r, based on the best-fi t 
values and the errors obtained b y ICalverlev et al.l (|201ll ) . 
In practice, ICalverlev et al.l l|201ll ) quote errors on the log- 
arithmic quasar luminosity In Lq. For simplicity, we mul- 
tiply the overall likelihood p{Rs , fr , x^^) , as defined in 
the previous section, by a Gaussian in In/r, centered at 
In /r = (0.982, 0.536, 1.58) for J1030-h0524, J1623-f 3112 and 
J1148-I-5251, respectively, and with the same uncertainty of 
fin/ ~ 0.31 for all three quasars. 

(5) MH[07 quoted the values of the parameters where 
the KS probabilities fell to l/3rd, l/9th, and l/27th of their 
peak values. For a ID Gaussian, these would correspond 
to ~ 1.5, « 2, and ~ 2.5a confidence limits. The proce- 
dures discussed above allow us to integrate the likelihoods 
within the fixed contours, and to compute actual confidence 
limits in three dimensions. This is a particularly important 
change: due to the long non-Gaussian tails, and having a 
three-dimensional parameter space, we find that the confi- 
dence levels differ significantly from the simple ID Gaussian 
expectation. 

In Table [TJ we show the parameter combinations where 
the likelihoods peak, as well as the best-fit models returned 
by the bootstrapping procedure. The first row for each 
quasar lists the results for the uniform-ionization case. For 
comparison, the values from MH07 are listed in the sec- 
ond row. The likelihoods peak at high neutral fractions for 
all three quasars, though none precisely at x^i^ = 1.0. A 
notable difference from MH07 can be seen in the case of 
J1148-I-5251: the KS probability peaks at a higher neutral 
fraction (a^Hi^ ~ 0.88 vs. 0.16), though this value is still con- 
sistent with the l/3rd probability contour of MH07. We also 
find significantly higher Rs values for both J 1030-1-0524 and 
J1623-|~3112 compared to MH07, caused by the increase in 
the intrinsic redshifts of these two quasars. This also leads to 
a factor of ~ 2 — 3 larger inferred quasar luminosities, which 
is not surprising: the material at the same observed wave- 
length is now farther from the quasar, so a larger intrinsic 
luminosity is needed to produce similar absorption. 

When we use the same Qks statistic as MH07, we 
further recover very similar lower limits on Xhi'^, at- 
tributable to GP-damping wings for both J1030-I-0524 and 
J1623-I-3112. As in MH07, we also find that low values of 
x^^j^ can not be excluded for J 11484-5251. 

We then switch to the new KS probability densities and 
integrate these over the parameters to obtain confidence lev- 
els. Our 68% and 95% two-dimensional confidence contours 
for each of the three quasars are shown in Figure (2] The con- 
straints shown in each panel have been marginalized over the 



third parameter. Focusing on the further marginalized con- 
straint on the single parameter Sfn'^, we obtain a clear lower 
fimit for J1148-H5251 and J1623-h3112 at about i^g'^ > 0.04 
while J 1030-1-0524 has a somewhat stronger lower- limit at 
xWj^ >Q.07 . 



3.2 Patchy Reionization 

We next fit the spectra using our patchy reionization mod- 
els. Unless stated otherwise, all constraints quoted below 
correspond to 95% confidence limits. 



3.2.1 Results with a Umform Prior on Rs 

In Figure [3l we present the joint 2D 68% and 95% confi- 
dence contours for each quasar. As in the uniform reioniza- 
tion models, the best fits are located at high neutral frac- 
tions (x^^j^ ^ 0.38). The most conspicuous change, however, 
is that the contours expand to significantly lower values of 
^m'^) with the lowest value ^hi'^ = 0.01 still falling inside 
the 95% CL on the — fr plane. The single-parameter 

marginalized constraint yields essentially the same lower 
limit, weakened to Xhj^ > 0.02 for all three quasars. 



3.2.2 Adding a Physical Rs Prior 



As described in § 12.11 we shifted the ionization field in the 
patchy-reionization simulation, in order to ensure that the 
quasar's ionization front can encounter a neutral pixel at the 
pre-specified location Rs. This was necessary to build up 
the statistics of the pixel optical depth PDF, and to be able 
to treat Rs as one of our free parameters. However, as de- 
scribed above, these shifts make fits to the observed quasar 
spectra unfair: clearly, it is less likely to find a neutral pixel 
within a fixed radius Rs in a highly ionized universe than in 
a more neutral universe. Thus, we next add in our analysis 
the prior probability for the LOS to intersect a neutral patch 
within Rs, taken from the patchy reionization simulations. 
These probabilities are shown in Figure [l] 

In Figure[5]we present the confidence regions with these 
prior probabilities included. Unsurprisingly, adding the pri- 
ors has the effect of excluding models with lower neutral 
fractions. For each of the three quasars, the 95% contour 
excludes models with Sjn^ < 0.05. The single-parameter 
marginalized constraint yields ~ 7-9 times tighter limits, 
with J1148+5251 and J1623-(-3112 at x}f;^ > 0.14 and 
J1030-|-0524at x^^^ > 0.11. We consider this to be the 
"fairest" statistic, and present these constraints as the main 
result of this work. 



4 DISCUSSION 

4.1 The Detection of the GP Damping Wing 

The main result of this paper is the strong lower limit on the 
mean neutral hydrogen fraction, Khj^ > 0.1. The mean GP 
damping wing absorption profiles, averaged over all LOSs 
in our best-fit models (in the patchy reionization case, with 
the Rs prior included, and with the best-fits correspond- 
ing to the peak of the P{Dks) likelihood) are shown, for 
illustrative purposes, by the dot-dashed (green) curves in 
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Figure [TrI We wish to confirm that our constraints arise 
from the statistical preference for these GP damping wings 
in the observed spectra, and not because of some other fea- 
ture0 To test this explicitly, we first re-fit the spectra, but 
this time ignoring the absorption by the GP damping wing. 
At the location of the best-fit model (i.e. with Rs and fr 
fixed at the values listed in Table[l|, the KS probabilities de- 
crease by more than an order of magnitude for J1030-I-0524 
and J1623-I-3112 and more than two orders of magnitude for 
J1148-I-5251. This demonstrates the sensitivity of the data 
to the presence/absence of the damping wings. 

Even if the data is sensitive to the damping wings, how- 
ever, our fits might not be sensitive to their presence, owing 
to degeneracies with other parameters. Indeed, we found 
that when we allow Rs and fr to fioat, the best-fit value 
of /r was lower for J1623+3112 and J1148-I-5251 than in 
our fiducial models which include a damping wing. This is 
unsurprising - when the quasar flux is lower, the resonant 
opacity is increased and can compensate for the loss of the 
damping-wing absorption. However, in these damping-wing- 
less best-flt models, the peak probabilities are still factors 
of ~6-ll below those in the original best-flt models. This 
demonstrates that our other parameters, e.g. a lower quasar 
ionizing flux, cannot fully mimic the presence of a GP damp- 
ing wing. 



4.2 Parametric Bootstrapped Confidence 

Contours and Scatter in the Damping Wing 

The confldence contours shown in Figures [2l [3] and [5] rely 
on the assumption of an underlying flat prior on log a^Hi'^ i^^ 
well as either a flat or a physical prior on Rs, and a Gaus- 
sian on log /r). Since the choice of the fiat prior on log 
is essentially arbitrary, we have replicated our entire analy- 
sis with an alternate, parametric bootstrapping procedure. 
Namely, we repeat the analysis discussed in § 12.4.21 another 
2 X 10* times, replacing the observed spectrum with the 
mock spectra generated along each of the 2 x 10'* LOSs in 
our best-fit model. This results in a set of 2 x 10* new best-fit 
parameter combinations, which, in general, differ from the 
input values. The number of the new best-fit values in our 
3D parameter space then directly yields the joint confidence 
levels. This method does not require an explicit prior on how 
the parameters are distribute d, and is a comm on technique 
to estimate confidence levels l|Press et al.lll992l V 

The results of this procedure are shown in Figure [6l 
The analysis includes the physical Rs prior from Figure [4] 
(here used only to identify the initial best-fit models from 
which the LOSs are bootstrapped), and so this figure is to be 
compared to Figure (5] As this comparison reveals, the boot- 
strapping returns noticeably tighter contours. The single- 

^ Interestingly, the damping wings produce some absorption, at 
the 5-10% level, extending to the red side of the quasar's hya 
emission line. This absorption was neglected in our modeling of 
the quasar's intrinsic spectrum. Taking this absorption into ac- 
count would increase the intrinsic flux; since an overall increase 
is degenerate with /r, we do not expect this to change our con- 
clusions on the inferred damping wing. 

For example, had we not blacked out our model spectra in the 
region R > Rs, we would have been sensitive to the absorption 
statistics in the IGM, as mentioned in a footnote in § 12.21 



parameter constraints on a;}^'" have tightened to xj^^ > 
0.48 for J1030-f0524 and ^hI^ > 0.62 for J1623-I-3112 
and J1148+5251. These constraints are weaker without the 
Rs prior, in which case the lower limit for J1030-I-0524 is 
^m^ > 0-11, J1623-H3112 is x^^^^ > 0.07, and J1148-H5251 
is > 0.14. The lower limits derived for all three quasars 

in the corresponding Bayesian approaches are much weaker. 

The large differences between the bootstrapped and the 
Bayesian confidence regions has a simple physical interpre- 
tation. The bootstrapping procedure samples the LOSs in 
our original best-flt models, and therefore the extent of the 
confldence regions reflect the LOS-to-LOS variations of the 
mock spectra in these best-flt models. For each quasar, this 
best-flt model corresponds to a nearly neutral IGM. Most 
importantly, in these nearly-neutral models, there is rela- 
tively little scatter in the strength of the damping wings - 
every LOS in each model contains a wide swath of neutral 
IGM. In contrast, low-x]^^ models contain only relatively 
small and isolated neutral islands. For example, in the model 
with xj^^ = 0.01, we find that neutral patches have a me- 
dian diameter of only ~ 7.5 Mpc, and are typic ally sepa- 
rated by distances of order ~ 100 Mpc (see also iMesingen 
|2010D . It is therefore very unlikely to find any LOS in our 
high-x^^ models, that is best fit by the narrow and weak 
damping wings in low-s]^^ models. This results in the tight 
constraints generated by the bootstrapping procedure. 

The situation is quite different when we estimate con- 
fidence levels directly from the Bayesian likelihoods. In this 
case, the probability density at each point in the parameter 
space is estimated using the LOS-to-LOS variations m that 
model, which generate the r CPDFs. As a result, low-x^^ 
models are ruled out at a much weaker significance: these 
models have a large scatter. For example, in the model with 
^m'^ — 0.01, we find that ~ 5% of the sightlines contain 
much larger neutral patches, with a diameter > 35 Mpc. We 
conclude, therefore, that the \ow-x^^ models have a rela- 
tively wide probability distribution of the KS-distance Dks, 
and are therefore difficult to rule out at a high significance. 

In principle, the scatter in resonant absorption may be 
important as well: the sightlines in the low-x^*^ models that 
look more similar to the QSO (i.e. with Dks < DKS.obs) 
might be those that have unusually high neutral density in- 
side the HII region. To study this, we compare a sub-sample 
of "good-fit" LOSs in our x]^^^ = 0.01 model (i.e. LOSs 
having low values of Dks < DKS.obs) to the full sample of 
2 X 10* sightlines. The Iow-Dks sub-sample has a notice- 
ably stronger damping wing than the full sample. In con- 
trast, the median resonant optical depths between the two 
samples are indistinguishable. We therefore confirm that the 
damping wing statistics are driving our conclusions and that 
the scatter in the damping wings make the Bayesian and the 
bootstrapping confidence limits different. These conclusions 
also show that the scatter in the damping wing strengths and 
shapes has a strong dependence on the average global ion- 
ization ^Hi'^, and imp ly that this can provide an ad ditional 
probe of reionization (|Mesinger fc FurlanettdboOSl 'l. 

Apart from the three quasars analyzed here (as well 
as the z = 7.1 source ULAS J1120-^0641) with fuU GP 
troughs, currently th e re are ^ 20 quasar s known at z > 6 
(jWillott et al.l l2010bl : ICarilli et all bOld l. The sample~is 
quite heterogeneous, and many quasars are too faint for use- 
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ful spectroscopy. However, there are a handful of sources 
with spectra whose quaUty is similar to the quasars we an- 
alyzed, and yet they do not show a fully black GP trough. 
This may place an upper limit on a;J§'^, since open sight- 
lines, with no neutral patches, become rare in a nearly neu- 
tral universe. 



4.3 Parameter Dependencies and Degeneracies 

The confidence contours presented in Figures [2] and [3] illus- 
trate a number of expected basic trends. Most importantly 
for our purposes, as the value of x^i^ increases, the GP 
damping wing becomes stronger, and, in particular, domi- 
nates further blueward. In the patchy models, the depen- 
dence on aSjji'^ arises from the width of the neutral slab(s) 
beyond R > Rs- This is different from the uniformly ion- 
ized models, in which the damping wing optical depth sim- 
ply scales linearly with a^Hi'^, at every frequency. As /r 
increases, the resonant absorption inside R < Rs is sup- 
pressed, and, as Rs increases, the transmission window gen- 
erally becomes larger. Degeneracies are therefore expected: a 
larger Rs can be partly compensated by a stronger damping 
wing (larger ijn^) and/or resonant absorption (a smaller 
/r). Likewise, although their wavelength-dependence is dif- 
ferent, in general, a weaker damping wing can be compen- 
sated, at least partly, by stronger resonant absorption. These 
trends are indeed visible in Figures [2] and [3] (although tem- 
pered by the priors on /r and on -Rs). 



4.4 Damped Lyman Alpha Systems 

We have thus far considered GP absorption from neutral 
material in the IGM (i.e. at densities near the mean back- 
ground density) . Here we consider an alternative scenario, in 
which the IGM is highly ionized, and the apparently black 
GP trough in the quasar spectrum is caused instead by a 
single discrete absorber, i.e. a DLA, along the LOS. 

We consider DLAs with the range of column densities 
of iVm = {10'^5x 10'^102^5x 102°,102\5x lO^'jcm-^, 
placed at the same distances Rs as we used for the edges of 
the surrounding HII regions in our fiducial analysis above. 
We create mock absorption spectra as outlined above, but 
replace the IGM damping wings from our reionization sim- 
ulations with DLA damping wings. There are no neutral 
patches in this model, and the ionization field is taken to 
be the sum of the quasar's ionizing fiux, parameterized by 
/r as before, and the fiux of the background galaxies, which 
is taken from the semi-numerical simulations as described 
above. Also as before, the fiux is conservatively assumed 
to be zero beyond the DLA's location {R > Rs)- This as- 
sumes that the DLA casts a full shadow of the quasar along 
the LOS, and makes our constraints independent of the ab- 
sorption statistics of the background IGM, away from the 
quasar. 

In Figure[71 we present joint two-dimensional confidence 
contours in this DLA model. As the figure shows, all three 
quasars require DLAs with very large column densities, with 
the same best-fit value of A^hi = 5 x 10^° cm^. As Table [1] 
further shows, the bootstrapping procedure yields 68% CL 
lower limits on A''hi (marginalized over both Rs and /r) 
that range between ^ 2 x 10^^ — 3 x 10^° cm^. Interestingly, 



Table [T] also shows that for J 1030+0524, the lower limits on 
A^Hi from the Bayesian and the bootstrapping analyses are 
similar, but for J1148+5251 and J1623+3112, the Bayesian 
limits are much weaker and closer to A^hi ~ 10^^ cm^^. This 
is similar to the differences we found for the lower limits 
on Xjn'^- However, unlike in the patchy reionization mod- 
els, the DLA damping wings have, by definition, no scatter 
to account for this. To investigate the origin of this differ- 
ence, we compared a sub-sample of "good-fit" LOSs — i.e. 
LOSs having low values of Dks < DKS,obs — in our low- 
A'hi model (with A^hi = 10^^ cm~^) to the full sample of 
2 X 10'' sightlines. The low-Dxs sub-sample had noticeably 
stronger resonant absorption than the full sample. We thus 
conclude that the difference between the Bayesian and boot- 
strapped lower limits for the DLA column densities is due 
to the scatter in the resonant absorption. This scatter is 
relatively more important in the low-Aini models, and can 
produce LOSs that mimic stronger DLAs. 

Most importantly, these results show that DLAs along 
the LOS can provide a plausible explanation of the spec- 
tra themselves. They could produce the black GP troughs, 
and, as shown in Table [1] the probabilities of the best-fit 
models are comparable to those in the models with incom- 
plete reionization. However, finding DLAs with the high re- 
quired column densities, relatively near the quasars, is im- 
probable. Because of their rarity, the abundance of DLAs 
at high redshift is poorly known. Nevertheless, a mod- 
est extrapolat ion of existing DLA mea s urements between 
4 < z < 6 l|Prochaska fc Wolfd l2009l : ISongaila fc Cowid 
I2OIOI ) to 2 « 6.5 imply the abundance ~ 0.05 DLAs with 
column densities Ahi ^ 2 x 10 ^° cm^ per unit redshif t 
(see, e.g.. Fig. 8 and Eq. 4 in ISongaila fc Cowid bOld ). 
In our case, ~ 50Mpc corresponds to Az ~ 0.1, im- 
plying that an abundance of ~ 10 DLAs per unit red- 
shift is needed: this is more than two orders of magnitude 
higher than the above extrapolation suggests. In principle, 
the abundance of DLAs around quasars could be atypi- 
cally high; however, at distances as far as 50 Mpc from 
the quasar, these biases should be negligible. Finally, the 
high-resolution spectra of these quasars show no evidence 
at the required dista nces of the low-ionization metal lines 
typical of DLAs (e.g. | Rvan- Weber. Pettini fc Madaul I2OO6I : 



iRvan- Weber et al.l l2009l : ISimcod I2OO6I : G. Becker, private 
communication 2010). 



4.5 The Impact of Ionizing Radiation from 
Galaxies 

We turn now to the question of whether the ionizing radia- 
tion from the galaxies in the quasar's ionized region is impor- 
tant for our conclusions. Generally, we expect that galaxies 
could be important for a combination of small /r (so that 
the quasar does not outshine the galaxies), small x^i^ or 
A^'hi (so that the damping wing absorption, which is inde- 
pendent of the galaxy background, does not dominate over 
the resonant absorption), and large Rs (since the galaxy 
fiux is only weakly dependent on R within R < Rs, while 
the quasar fiux drops as ~ ^/R^) 

To investigate the importance of galaxies more explic- 
itly, we have re-run our analysis, with the ionizing radiation 
from the galaxies artificially turned off inside R < Rs- We 
find that for larger values of fr the galaxy flux can be ig- 
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nored without much consequence. Furthermore, the presence 
of the galaxy flux has almost no effect on the best-flt loca- 
tions or the peak likelihoods, with the impact restricted to 
the "periphery" of our parameter space, as expected above. 
In general, the contours tend to expand toward larger Rs 
and smaller aij^^ when the galaxy flux is not modeled. 

As an example, we find that DLA models for 
J1623-I-3112 are most affected by the galaxy flux. The prob- 
ability density at (iVni = lO^^cm"^, Rs = 47 Mpc), inte- 
grated over /r, is ~ 3 times higher when the galaxy con- 
tribution to the ionization radiation is neglected. The ef- 
fect is reversed in (high-A^ni, low-J?s) models where, for 
example, in the same quasar for A^'hi ~ lO^^cm"^ and 
Rs = 37 Mpc, the probability density (again integrated 
over /r) is ~ 10 times lower without galaxy fluxes. Nev- 
ertheless, the marginalized ID constraints remain robust: 
without galaxies, following the format of Table [T] the con- 
straints are (/p, i?s, log iVm) = (1.61^;°;^, 39lf , 34, 20.0155;^^) 
with max (p) = 0.21. 

4.6 The Impact of Self-Shielded Absorption 
Systems along the LOS 

In this section, we consider the possibility that attenuation 
along the LOS is substantial, causing the QSO's intensity 
to fall off more rapidly with distance than the approximate 
oc profile in equation ((3]). As mentioned above, a DLA 
or strong LLS along the sightline would be rare, and likely 
detectable in HIRES spectra. However, several lower column 
density systems aligning along the LOS could in principle 
strongly attenuate the fiux, thereby mimicking a damping 
wing signature. 

How degenerate are these two signatures? To answer 
this question, we perform an additional set of runs, repeat- 
ing our fiducial analysis above, except without a damping 
wing, but instead attenuating the QSO fiux by a factor of 
exp(— r/Ag). We vary the effective LOS mean free path, 
Aq, from 60 Mpc down to 5 Mpc. We find that in all of 
these runs, the peak probabilities are reduced from our fidu- 
cial models with damping wings (Table [T]) . Specifically, for 
J1623-h3112 / J1030-h0524 / J1148-h5251, the best-flt mod- 
els have Xq = 50 / 40 / 15 Mpc, with max (p) lower than the 
peak probabilities in the flducial, damping-wing-included re- 
sults by 76 / 97 / 34 %. 

To further explore the degeneracy between an altered 
(i.e. steeper) resonant absorption proflle and our fldu- 
cial damping-wing absorption, we construct an unphysical, 
"worst-case" model. We again have no damping wings, but 
instead we adjust the fiux, in each pixel, such that the aver- 
age resonant optical depth precisely matches the total (res- 
onant+damping) optical depth in our best-fit models from 
Table [1] We thereby allow an entirely arbitrary and con- 
trived profile of the flux along the line-of-sight. (Indeed, in 
some cases, this proflle is unphysical, as it rises with dis- 
tance away from the quasar.) Even in this worst-case sce- 
nario, we find that the peak probabilities are still lower 
than in our fiducial cases, by 23% for J1148-I-5251, 11% for 
J1030-h0524, and 22% for J1623-H3112. This preference for 
the presence of a damping wing implies that the spectral 
fits are (mildly) sensitive to the smoothness of the damping 
wing profile (compared with the resonant absorption profile, 
which has significant pixel-to-pixel variations). 



From both of the above tests, we conclude that even if 
there is unaccounted-for absorption of the QSO flux along 
the LOS, the observed spectra still prefer an IGM damping 
wing profile, at the tens of percents level. 



5 CONCLUSIONS 

The main result of this paper is the lower limit on the mean 
neutral hydrogen fraction, Spji^ ^^-^ (^^ 95% CL), inferred 
from the spectra of three z > Q quasars. Lower limits were 
previously obtained for the same three quasars by MH07, us- 
ing simpler assumptions, and applying a less sophisticated 
statistical analysis. Our new results strengthen these pre- 
vious limits. In particular, we have found that low neutral 
fractions are ruled out after relaxing the major assumption 
in MH07, of a uniform ionizing background. Unlike the lower 
limit obtained by MH07, which was driven by the pixel op- 
tical depth statistics in the spectral fitting alone, our new 
limit is driven by our additional modeling of the spatial dis- 
tribution of neutral patches. We note that if we conserva- 
tively allow the size of the surrounding HII region to be a 
fully free parameter, independent of aSfn^, our constraints 
weaken to s^'f > 0.03. 

Our results arise through a statistical preference for a 
GP damping wing absorption component in the spectra. We 
confirm that peak likelihoods drop by a factor of several to 
more than an order of magnitude when the spectra are mod- 
eled without the damping wing component, implying the 
preference for significant neutral hydrogen in the IGM. Al- 
though such a damping wing could arise from a high-column 
density (TVhi > few x lO^^cm"^) DLA, such systems are very 
rare; furthermore, high-resolution spectra show no evidence, 
at the required distances, of the low-ionization metal lines 
typical of DLAs. In principle, the imprint of the damping 
wing could be mimicked by an alternate (not included in 
our fiducial model) evolution of the resonant absorption with 
distance from the quasar. We have shown, using a contrived 
worse-case toy-model, that the statistical preference for a 
damping wing can indeed decrease, but that it can not fully 
go away, even in this contrived model. Our analysis implies 
that reionization is not yet complete at 2: ~ 6.2. 

Finally, we find different results when a Bayesian or 
a parametric bootstrapping method is used to estimate the 
confidence levels on our model parameters. We attribute this 
difference to the fact that the two methods are sensitive to 
the sightline-to-sightline scatter in the absorption statistics 
in the low- and the high-a;^'^ models, respectively. The 
scatter, in particular, in the GP damping wing strengths 
and shapes has a strong dependence on the average global 
ionization - as the global neutral fraction increases, the scat- 
ter is diminished. Our results demonstrate that this scatter 
could provide a useful additional diagnostic of the ioniza- 
tion topology of the IGM, when applied simultaneously to 
a larger quasar sample in the future. 
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Table 1. The table shows the best-fit parameters in different models for the three quasar spectra (along with their 95% errors when 
appropriate, marginalized over the other two parameters in each case). For each quasar, the five rows describe the following models: (i) 
uniform reionization (fj I3.1|l . (ii) the earlier results of MH07, for reference, (iii) patchy reionization fS I3.2.1ll . (iv) patchy reionization with 
a physical prior added for the distance Rg to the nearest neutral patch (fj 13. 2.21 1. and (v) a model in which the GP trough is caused 
by a DLA (fj I4.4| l. For each model, the four rows show (i) the best-fit parameters that maximize the probability based on the KS test 
(the values listed are for (/r, iJg/Mpc, x^'^) in the reionization models, and (/r, i?s/Mpc, log Nm/cm"^) in the DLA models), (ii) 
the value of this peak probability p; (iii) the best-fit parameters from the bootstrapping procedure (except MH07 did not perform any 
bootstrapping), and (iv) the peak number (Ni,f) of bootstrap best-fits out of 2 X 10** trials in this model. *Note that we report 68% 
confidence intervals for logA^Hi. 
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Figure 1. The continuum— normalized spectra of the three SDSS 
quasars used in this paper. The dotted (blue) vertical lines in 
each panel mark the boundaries of the three wide wavelengths 
bins used in the analysis, while the dashed (red) lines show the 
wavelengths corresponding to various assumed sizes (Rs) of the 
quasar's ionized region explored in our models. The dot-dashed 
(green) curves indicate, for illustrative purposes only, the trans- 
mission attributed to the GP damping wing alone. These are the 
mean damping wings in the best-fit (without an Rg prior) models 
for each quasar, but our actual analysis includes the LOS-to-LOS 
scatter in both the damping wing and the resonant absorption). 
The magenta arrow on the right of each panel marks the location 
of the hya line at each quasar's redshift. 
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Figure 2. Confidence regions under tiie assumption that tiie IGM 
is ionized by a spatially uniform background. The three rows show 
the results for the three different quasars. In each row, the three 
different panels show the confidence contours in three different 2D 
planes, marginalized over the third parameter. The contours en- 
close 68 and 95% of the marginalized 2D likelihood. The location 
of the best-fit model is marked by an "x" in each panel. 
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Figure 3. Confidence regions as in Figure[2l except in our patchy 
reionization models. 
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Figure 4. Fraction of the simulated lines-of-sight (plos) that 
intersect a neutral pixel within a given distance away from the 
quasar's host halo, as labeled. These fractions are obtained from 
our patchy reionization models, and shown as a function of the 
global average neutral fraction S^'^. The figure illustrates that 
intersecting a neutral patch within 35-60 (comoving) Mpc from 
the quasar is much less likely if the IGM is highly ionized. 
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Figure 5. Confidence regions as in Figure [S] except priors were 
added for the probability of finding a neutral patch at Rs , taken 
from Figure |4] 
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J1148 




Figure 6. Confidence regions as in FigurefS] except that the para- 
metric bootstrapping procedure, discussed in S I2.4.2I was used to 
estimate the confidence contours. 
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